/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package provider.impl.genomemap.lc;

import commons.core.LikelihoodComputer;
import genomemap.data.TAB08Data;

/**
 * @author Susanta Tewari
 * @since Aug 21, 2011
 */
final class TAB08LCRLA extends TAB08LCSFA implements LikelihoodComputer<TAB08Data> {

    /**
     * allele pattern for the tetrad at each locus
     */
    private final byte[][] tetrad_pattern = {

        { 1, 1, 0, 0 }, { 0, 1, 1, 0 }, { 1, 0, 1, 0 }, { 1, 0, 0, 1 }, { 0, 1, 0, 1 },
        { 0, 0, 1, 1 }
    };

    /**
     * matching pattern for the genotype with the tetrad for each previous interval in the iteration
     * of the recursive linking algorithm
     */
    private final byte[][] match_pattern = {

        { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { 0, 0, 0, 1 },
        { 1, 1, 0, 0 }, { 1, 0, 1, 0 }, { 1, 0, 0, 1 }, { 0, 1, 1, 0 }, { 0, 1, 0, 1 },
        { 0, 0, 1, 1 }
    };

    @Override
    double[] implementHMM(final byte[][] f, final double[] c) throws InterruptedException {

        final int loci = c.length + 1;

        if (loci == 2) {
            return super.implementHMM(f, c);
        }

        final double[] lhood = new double[f.length];

        /* a general lhood for 11 matching pattern and 6 tetrad pattern */
        double[][][] old_general_lhood = null;

        /* interval position not index. It runs from loci-1 to 1 */
        for (int intervalPos = loci - 1; intervalPos > 0; intervalPos--) {

            final double[][][] general_lhood = new double[f.length][11][6];

            for (int i = 0; i < match_pattern.length; i++) {

                for (int j = 0; j < tetrad_pattern.length; j++) {

                    final ExchangeConfig x = new ExchangeConfig();

                    do {

                        for (int index = 0; index < f.length; index++) {

                            final byte[] strand = tetrad_pattern[j].clone();

                            modifyStrand(strand, x.getI());
                            modifyStrand(strand, x.getJ());

                            /* match genotype (f) with this tetrad for the interval */
                            final byte[] match = new byte[4];

                            for (int k = 0; k < 4; k++) {

                                if ((tetrad_pattern[j][k] == f[index][intervalPos - 1])
                                        && (strand[k] == f[index][intervalPos])) {
                                    match[k] = 1;
                                }
                            }

                            /* update match with the previous interval tetrad matching pattern */
                            if (intervalPos != 1) {

                                for (int k = 0; k < match.length; k++) {
                                    match[k] *= match_pattern[i][k];
                                }
                            }

                            /* if last interval */
                            if (intervalPos == loci - 1) {

                                /* sum the matches */
                                byte sum = 0;

                                for (final byte val : match) {
                                    sum += val;
                                }

                                final double marginal_prob =
                                    getMarginalProb(new double[] { c[intervalPos - 1] },
                                                    new ExchangeConfig[] { x });

                                general_lhood[index][i][j] += sum / 4.0 * marginal_prob;
                            } else {

                                /* look up k for 'match' in old_general_lhood */
                                final int i_lkup = indexOf(match_pattern, match);
                                final int j_lkup = indexOf(tetrad_pattern, strand);

                                if (intervalPos != 1) {

                                    general_lhood[index][i][j] +=
                                        old_general_lhood[index][i_lkup][j_lkup]
                                        * getMarginalProb(new double[] { c[intervalPos - 1] },
                                                          new ExchangeConfig[] { x });

                                } else {

                                    lhood[index] +=
                                        old_general_lhood[index][i_lkup][j_lkup]
                                        * getMarginalProb(new double[] { c[intervalPos - 1] },
                                                          new ExchangeConfig[] { x });
                                }
                            }
                        }

                    } while (getCrossoverPattern(x));

                    if (intervalPos == 1) {
                        return lhood;
                    }
                }
            }

            old_general_lhood = general_lhood.clone();
        }

        throw new Error("control should not come here.");
    }

    /**
     * Method description
     *
     *
     * @param array1
     * @param array2
     *
     * @return
     */
    private int indexOf(final byte[][] array1, final byte[] array2) {

        boolean flag = true;

        for (int i = 0; i < array1.length; i++) {

            flag = true;

            for (int j = 0; j < array1[i].length; j++) {

                if (array1[i][j] != array2[j]) {

                    flag = false;

                    break;
                }
            }

            if (flag) {
                return i;
            }
        }

        return -1;
    }
}
